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NEW EPITAXIAL THIN FILM MODELS AND NUMERICAL APPROXIMATION 


WENBIN CHEN*, ZHENHUA CHENt, JIN CHENG*, AND YANQIU WANG® 


Abstract. This paper concerns new continuum phenomenological model for epitaxial thin film growth with three different 
forms of the Ehrlich-Schwoebel current. Two of these forms were first proposed by Politi and Villain m and then studied by 
Evans, Thiel and Bartelt [5]. The other one is completely new. Following the techniques used in Li and Liu m, we present 
rigorous analysis of the well-posedness, regularity and time stability for the new model. We also studied both the global and the 
local behavior of the surface roughness in the growth process. The new model differs from other known models in that it features a 
linear convex part and a nonlinear concave part, and thus by using a convex-concave time splitting scheme, one can naturally build 
unconditionally stable semi-implicit numerical discretizations with linear implicit parts, which is much easier to implement than 
conventional models requiring nonlinear implicit parts. Despite this fundamental difference in the model, numerical experiments 
show that the nonlinear morphological instability of the new model agrees well with results of other models published in 
which indicates that the new model correctly captures the essential morphological states in the thin film growth process. 


Key words, epitaxial thin film growth, Ehrlich-Schwoebel effect, convex-concave splitting method, semi-implicit time 
discretization. 
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1. Introduction. In epitaxial thin film growth, the phenomenological continuum evolution of film height 
h{x,t) at lateral position x G C and time t G (0,r) is governed by the equation [T^ 

dth = C-w-J, (1.1) 

where is a given function related to the deposition rate, and J(x, t) is the lateral mass current of 

adatoms across the film surface. The current J consists of an equilibrium (EQ) part and a non-equilibrium 
(NE) part, denoted by J = Jbq + Jab- For the equilibrium part, we adopt the linearized model of Mullins 
El and set Jeq = KEQ'^{Ah), where the constant Keq > 0 is usually very small. The more interesting 
non-equilibrium surface current Jne depicts the interaction of adatoms with surface steps, and here we follow 
the model presented by Evans, Thiel and Bartelt [5]: 

Jnb = Jdf + Jbs + ^ relax, 

where J_df = ~ 7 V/i, with constant 7 > 0, is the stabilizing downward tunneling (DF) current; 3es is the 
de-stabilizing uphill Ehrlich-Schwoebel (ES) current, which will be discussed in further details later; and 
^RELAX = FV(Ah), with constant k > 0, is a phenomenological relaxation current artificially added when 
Jeq ~ 0, in order to counteract the increasingly violent unstable behavior caused by Jfs- Mathematically, 
one can combine 3relax with 3eq to get 3eq + 3relax = e^V(Ah), with = k -I- Keq- 

Now let us examine 3es, which models the Ehrlich-Schwoebel (ES) effect. The ES effect states that 
adatoms must overcome a higher energy barrier in order to attach to a step from an upper terrace than from 
a lower terrace. Thus it favors an uphill current and consequently causes the formation as well as steepening 
of mounds OiniEolEl]. Due to its nonlinear nature, the ES current brings interesting surface morphological 
instability, but imposes difficulty upon the mathematical analysis. To our knowledge, there exist three ES 
models which have been mathematically investigated in terms of well-posedness and properties of the solution: 

1. Infinite ES barrier model proposed in [53] with ES current 3i es = \^h \2 i 

2. Finite ES barrier model proposed in [9] with ES current Je,es = 

3. Finite ES barrier with slope selection model (see [T3]) with ES current Jess.es = (1 — |V/ip)V/i; 
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where | • | stands for the Euclidean length of a vector. Note that 3pss.ES and Jj.es can be viewed as 
approximations to 3f,es when |V/i| <C 1 and |V/i| ^ 1, respectively. In [T3], well-posedness and long time 
stability have been established for the two finite ES barrier models. We point out that a main difference, 
from the mathematical point of view, between the infinite ES barrier model and the finite ES barrier models 
is that, 3i.es is not continuous at V/i = 0, while Jf.bs and 3fss,es remain continuous for all V/i. This is 
why rigorous mathematical analysis in |13| only works for the two finite ES barrier models. 

In this paper, we are interested in a different form of ES current first proposed by Politi and Villain [18) . 
and later studied by Evans, Thiel and Bartelt [5]: 

f Vh Vh \ 

3i,es = cni — |„,| - , 1^,1 I , ( 1 - 2 ) 

Vp+|V/i| q+lVhfJ 

where oi > 0 and 0 < p < q are given parameters. There are more physical parameters involved in the 
real models in [^, and we only try to describe one simple but essential model here. Physical meanings and 
practical choices of these parameters will be presented in Section For now we focus on the mathematical 
properties of the ES current. 

When |V/i| ^ p, by eliminating high order terms one has 


JpBS = ai- 


P 


pq + {p + 9 )|V/i| + |V/ip 
Thus we introduce a variation of the ES current 

32,ES = OL2 


\Ih : 


ai- 


q-p 


Vh 


P J g +(p + 9 )|V/i|/p' 


V/i 


9+ (p + 9)|V/i|/p’ 


(1.3) 


with 02 = To our knowledge, this model of the ES current is brand new. 

Similarly, when |V/i| p, one has 

j ^ _ g-P __ f _ N_V/i_ 

_|_ 2|V/i|) + (g - p)|V/i| + |Vhp ^ ^ (g - p)|Vh| + |V/ip ’ 


because p(g + 2|V/i|) = p^ + p(g — p) + 2p|V/i| (g — p)|V/i| + |V/ip. This allows us to introduce another 
variation of the ES current 


Js.FS — 0^3 


V/i 

(g-p)|V/i| + |V/iP’ 


(1.4) 


with 03 = ai(g —p). A much simpler one-dimensional form of has been proposed and studied in [T51I5]. 

We believe this is the first time that the multi-dimensional form of is presented. 


The main purpose of this paper is to analyze mathematically the epitaxial thin film growth equation (1.1) 
with ES currents Jfc,FS) for A: = 1,2 and 3. Note that Js^fs is not continuous at V/i = 0 , while Ji,fs and 
J 2 ,FS are continuous for all V/i. In this sense, one may compare Js.fs with the infinite ES current J/,fs- 
Similarly, Ji_fs is comparable to the finite ES current without slope selection Jf,fS: and J 2 ,fs is comparable 
to the finite ES current with slope selection 3fss,es- Later it shall become clear that the models Jfe.FS, 
for k = 1,2,3, have built-in and significant differences from Jf,fs, Jfss.fs, and J/,fs in the mathematical 
analysis. Though interestingly, numerical results will show that they give very similar nonlinear morphological 
evolution processes, which is a good sign as they all model the same physical phenomenon. 

For simplicity, let fl be a rectangular domain and set the fl-periodic boundary condition on h. Following 
the previous discussions. Equation (O equipped with initial and boundary conditions can be written as 


dth = C + 7 A /1 — e^A'^h — V • JfSj in fl x (0, T], 

h{-,t) is ri-periodic for all t S [0,r], (1.5) 

h(x, 0 ) = ho(x) for all x € fi, 


where Jfs is chosen from Jfc.FS for A: = 1,2, 3. Here and throughout the rest of the paper, we shall only use 
subscript k when individual features of the ES current from different models are needed. Otherwise, the ES 
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current will simply be denoted as J es^ which can be any applicable choice from 3k ,es for k = 1,2, 3. Because 
of the singularity of at Vft. = 0, rigorous mathematical analysis in the rest of this paper will only be 

performed for Ji,_bs and 32,es- Though we still keep as an alternative option whenever applicable. For 

compatibility purpose, obviously ho and C should also be fl-periodic. 

Next, we introduce a surface roughness indicator and also argue that it suffices to study (1.5) under the 
assumption that C and ho are mean value free on 12. For simplicity, denote j^fdx= jiijy f^fdx for any given 
function /. Define the average height function h{t), for t G [0,T] hy h = hdx. Then, the surface roughness 
is estimated by m 


uj{t) = \ T \h{x, t) — dx, for t G [0, T], 


Similarly, denote C = f^C dx. By taking the average integral of the differential equation in (1.5) and using 


the ri-periodic boundary condition of h, it is clear that h satisfies an ordinary differential equation 


dth = (. 


Subtracting this equation from (1.5) and noticing that any spatial derivative of h is 0, one immediately gets 
dt{h -h) = {(-()+ 7 A(/i -h)- e^A2(/i -h)-W ■ 3Es{^{h - h))- 


In other words, h — h satisfies Equation (1.5), with C, in the right-hand side replaced by C ~ Ci the mean-value 
free component of C,. Thus studying h — h is equivalent to studying h with the assumption that C and ho are 
mean value free. In this case the surface roughness indicator becomes 

uj{t) = ^ j'Jh{K,t)\^ dx = -^||/i(-,t)||L2(n), forte [o,r]. 

Throughout the rest of this paper, we shall assume that and ho are mean value free on 12, and consequently 
so is h. Note that a typical example is C = 0. 

Using techniques similar to those in m, i e., the well-known Galerkin approximation and compactness 
argument approach of Lions m , we will establish the existence, uniqueness, and regularity of the weak solution 
to (1.5). The theoretical proof, although standard, relies heavily on particular properties of the ES current 


Jes- One of the main contribution of this paper is to establish these properties for 3k,ES-i with k = 1,2, and 
part of the properties for 3o,es- 

We will also establish global and local bounds for the surface roughness oj{t). The epitaxial thin film 
growth is in general a coarsening process, for which uj(t) is an important indicator. In the early stage of 
the growth, a typical rough-smooth-rough pattern mm is often observed. Hence theoretical and numerical 
study of uj{t) is important to the understanding of the surface morphological evolution. Besides the roughness 
indicator, the growth process is always energy driven in the sense that the dynamics is the gradient flow of a 
certain energy functional n 0 uni US [H]. We will show that the energy functional remains non-increasing 
with our ES current models, when the deposition rate (' = 0. 

Numerical discretization will be done using the convex-concave splitting technique. This technique was 
first proposed by Eyre to solve the Cahn-Hilliard and Allen-Cahn equations 0. Its main idea is to treat the 
“convex” part implicitly and the “concave” part explicitly in the time discretization. From another point of 
view, this is equivalent to solving a minimization problem of a strictly convex and coercive functional known 
as the modified energy functional [miMiiH]. Eyre’s convex-concave splitting scheme is first-order accurate in 
time and unconditionally stable. Later, higher order time schemes have been constructed using the similar 
idea um- For thin film epitaxial growth with ES currents 3e,es and Jess.es, the convex-concave splitting 
inevitably generates a nonlinear convex part [SHE], though alternative schemes with linear explicit parts can 
derived using other techniques mm- A significant advantage of the new ES models 3k,ESj for k = 1,2,3, 
is that, they naturally generate linear convex parts and nonlinear concave parts in the splitting. Hence the 
direct application of the convex-concave splitting technique will result in a linear problem to solve at each 
time step. Spatial discretization is done by a Fourier spectral Galerkin method. 
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The rest of the article is organized as follows. In Section we establish the existence, uniqueness, and 
regularity results of the weak solution to the model problem. In addition, bounds of the roughness indicator 
io{t) and analysis of the energy functional will also be given in this section. In Section a semi-implicit 
fully-discrete numerical schemes using the convex-concave splitting technique is presented. We show that the 
scheme is unconditionally stable. Convergence rate is also proved. In Section]^ we present numerical results 
which show similar morphological instability as results given in [18) . 


2. Well-posedness of the model problem. In this section, we study the well-posedness of Equation 
(1.5). As mentioned earlier, rigorous analysis will only be done when the ES current is taken to be either 
or 32^es- We shall first prove a few properties of the ES current in Section 2.1 in which we conveniently 
use a subscript fc = 1,2 to denote whether is taken to be Ji,_es or J 2 ,es, as the proof depends on the 
individual definitions of 3k,ES- K is worth to point out that Js.ss also possesses some similar properties, 
especially the most important convex-concave splitting one. This is why we do not want to completely leave it 
out, and the properties of will be mentioned in a separate remark. After these properties are established 
by fc-specific proofs, for simplicity we will drop the subscript k when the analysis does not depend on k. 


2.1. Properties of the function Jes- We start from 3k,ES for fc = 1,2. Note that 3k,ES depends 
solely on Vft.. It is convenient to view them as functions Jfc^Es(m) taking values at m = V/i. Moreover, by 
definition, we can write where 


$i(s) = Oi 


q-p 


{p + s){q + s)' 


and $ 2 ( 5 ) = 02 


q+ {p + q)s/p’ 


for all s > 0. 

Lemma 2.1. For all s > 0 and fc = 1,2, one has 


0 < $fc(s) <C, -C< < 0, 


where C is a positive general constant depending only on Uk, p, and q. 

Proof. The bounds for $fc(s), k = 1,2 are obvious, and the bounds for follows immediately from 


f \ 2(g-p)((p-hg)/2-hs) 

= -«i- , _^ 9 /_ , _N 9 - > 


2 (g -p) 


^ 2 ( 5 ) = - 0:2 


(p -I- s)2(g -I- s)2 

{v + q)/p 


(p + sY{q + s) 


> —Oi 


‘^■{q-p) 


p^q 


{q+{p + q)s/pY 


> — 02 - 


pq 


□ 

Another important observation is that both for k = 1,2, are gradient fields. Indeed, define 

functions Gk ■ —)■ K by 

— Gi(m) =pln(p-h |m|)-gln(g-h |m|) and —G 2 (m) =- ^|m|-hIn -h |m|V 

oi 02 p + q \p + qr \p + q J 

for all m G Now we examine the derivatives of Gfc(m) with respect to variable m. In order to distinguish 
such derivatives with the spatial derivatives, we use VEGfc(m) and V|;.Gfe(m) to denote the gradient and the 
Hessian of Gfc(m) with respect to m, while reserving the notation V and for gradient and Hessian with 
respect to the spatial variable x. 

Lemma 2.2. For k = 1,2, one has Gk G G^(K^). Their gradients satisfy 


VFGfc(m) = -4>fc(|m|)m, 


and their Hessians satisfy 

V|Gi(m) = -Oi 
V|.G2(m) = -02 


q-p 


[p + \ni\){q + \: 

P 

pq + {p + q)\ni 


-I + Oi 


1 


I + 02 


{p+ |m|)2 

p{p + q) 


(9+ |m|)^ 


m 0 m 


(pq + {p + q)\m\)‘^ 


m 


( 2 . 1 ) 
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where I is the 2x2 identity matrix and m(8)mzsa2x2 matrix defined by mm^, in which m is considered 
as a column vector. 

Proof. The proof is elementary. One only needs to use the fact that Vi?|m| = m/|m| and V p = 

S+Tinl“ (e+|L|)^ c> 0, to compute Vi;’Gfc(m) and V|Gfc(m). □ 

Corollary 1. We clearly have 3k, es € G^(K^) for fc = 1,2, and hence they are locally Lipschitz. 


Next we shall discuss the convex splitting of functions Gk{-), for A: = 1,2. We say a function is convex if its 
Hessian matrix is positive semi-definite everywhere, and concave if its Hessian matrix is negative semi-definite 
everywhere. It is not hard to see that 


Lemma 2.3. For all xi 7 ^ andx2^0i the function Gk{rn.) — ^Xk\'>^\'^; for k = 1,2, is concave. 

Proof. Note that the two eigenvalues of matrix m 0 m are 0 and |mp. By (2.1), it is clear that the two 
eigenvalues of VpGi(m) are 


Ai — —ai 
A 2 = —OLl 


q-p 


(P+ |m|)(g-)- |m|) 
q-p 


< 0 , 


+ ai 


1 


1 


|m| < ai 


1 


1 


= ai 

< 2ai 


(p-f |m|)(g-)-|m|) V(p-t|m|)2 (q + \ui\)^ J V(P+|m|)2 {q + \ni 

{q-p){p + q + 2\m\)\m\ f q-p \ f {p + q)/2 + \m\\^ ^ |m| 

p + |m 


m 


(p-f |m|)2(g-p |m|)2 
q-p 


= 2q!i 


(P+ |m|)(q-f |m| 


9+ |m| 


pq 


< Xi- 


This, combined with the fact that = Xi^: implies that Gi(m) — fyilmp is concave. Similarly, 

the two eigenvalues of V|-G 2 (m) are 


Ai — —02 
A2 = —02 


P 


pq+ (p + q)lnil 
P 


< 0 , 


02 “ 


p(p -f q) 


pq+ (p + q)lnil (pg-P (p-t g)|m|)2 


m = -02 


p 2 q 


(pg-f (p-f g)|m|)2 


< 0 . 


Hence G 2 (m) — Ax 2 |Lnp is concave. This completes the proof of the lemma. □ 

Corollary 2. The functions Gk{-), for k = 1,2, have the convex-concave splitting Gk{-) = Gk,+ {-) + 
Gk,-(-), where the convex and the concave parts are defined, respectively, by 

Gk,+ {m) = ^Xfclmp, Gfe _(m) = Gfc(m) - ^Xfolmp, 


for all xi > 2oi^ and X 2 > 0. 

The convex splitting and its properties are essential in theoretical analysis and the constructing of numer¬ 
ical schemes. In [13 [H], several bounds of the convex splitting for ES currents J f,es and J fss,es have been 
proved. Next, we shall prove similar bounds for ES currents with fc = 1, 2. 

Lemma 2.4. For any m G and k = 1,2, we have 

|Gfe(m)| < G(1 + |m|) and |V_FGfc(m)| < G, 

where G is a general constant depending only on otk, p and q. Moreover, all eigenvalues 0 /V|.Gfe(m) have 
absolute values bounded by G. In other words, the matrix 2-norm 0 /V|iGfc(m), denoted by |V|.Gfe(m)|, has 
bound 


|V|Gfc(m)| <G. 


Proof. The proof for |Gfc(m)| < G(1 -b |m|) follows immediately from the fact that ln(l -b s) < s for s > 0, 
while the proof of |VFGfc(m)| < G is elementary by Lemma 2.2 and the definition of 4’fc(-). Finally, the claim 
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about eigenvalues of V|.G/c(m) follows from the computation of these eigenvalues in the proof of Lemma 

□ 

Corollary 3. By Lemma\2.4\ and Corollary^ one immediately have 


2.3 


|Gfe,+ (m)| + |Gfc _(m)| < G(1 + |mp), 
|VFGfe,+ (m)| + iVj’Gfc _(m)| < G(1 + |m|), 
|V?.Gfe,+ (m)| + |V|Gfc,_(m)|<G, 


for any m S and k = 1,2. 

In addition, we also have the following lemma: 

Lemma 2.5. For any constant j3 > 0, there exists a > 0 such that 

Gfc(m) > —/?|m|^ — G/ 3 , for all m G and k = 1,2. 


Proof. By Lemma 2.4 and the Young’s inequality, one has 


Gfc(m) > -G(|m| + 1) > -G|m| - G > -/3|m|2 - ^ - G, 
where G is a positive constant. This completes the proof of the lemma. □ 


We have so far stated all properties of 3k,ES needed in the analysis of Equation (1.5). Note that these 
properties hold for both k = 1 and k = 2. It turns out that although not continuous at V/i = 0, also 

satisfy some of these properties. We summarize it in the following remark: 

Remark 4. Similar analysis shows that the same properties as presented in this subsection hold for 
as long as V/i stays away from 0. Below are the details. Define 


$3(s) = aa- 


1 


and —Ga(m) = — ln(g — p- 
ots 


{q-p)s 

Then, one has $a G GHK+), Ga G G(M2) n G2(M2 \{o}), and for all m G 
V_FGa(m) = -$a(|m|)m, 

1 , q — p + 2|m| 


V^Ga(m) = -aa- 


m (H) m 


Q!a- 


'(g-p)|m| + Imp ■ '-{{q-p)\m\ 

Moreover, Ga has the convex-concave splitting Ga = Ga,+ + Ga,- where 


Imp) 


G'3.+ (m) = ^Xalml' 


G3,_(m) = G3(m) - -X3|m|‘ 


for all xa > {q-%)'^ ■ When s or |m| stays away from 0, J 3 ,_es and Ga have similar bounds as in Lemmas 
2.4: 2.4 and Corollary^ but not when s —>■ 0 or 


m 


0 . 


2.1 


2.2. Weak solution to Equation (1.5). Due to the unboundedness of J 3 ,_es and G 3 mentioned in 
Remark]^ the analysis from here to the end of Sectiononly works for Jk,ES, with fc = 1,2. Using lemmas 
and corollaries proved in Section [2.1[ we no longer need to distinguish between k = 1 and fc = 2 in the analysis 
to be given. Therefore the subscript fc will be dropped for simplicity, i.e., without special mentioning, J_es, 
<!)(•), X and G(-) will be used with definitions taken to be either for fc = 1 or fc = 2 . Also, the convex splitting 
of G(-) defined in Corollary will simply be denoted by G+(-) and G_(-). Occasionally, the case fc = 3 will 
be discussed individually in remarks. 

In this subsection, we define what is a weak solution to Equation (1.5) and establish the existence, 
uniqueness as well as the regularity results of the weak solution. The analysis follows exactly the same 
framework presented in |13j . i.e.. Lions method m of first constructing a semi-discrete Galerkin spectral 
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approximation and then proving its convergence using a compactness argument, as this is currently the most 
efficient approach for the given problem. However, due to the different properties of 3ES^ there are still many 
essential differences between our analysis and the one in |13j . mainly in the proof of some inequalities. Thus 
we still present the entire proof for completeness, although readers may find the majority of notation and 
analysis are just borrowed from |13) . 


We first introduce the weak formulation of (1.5). Denote by for m > 0 and 1 < r < oo the 

H-periodic Sobolev space with indices m and r. When m = 0 and r < oo, the space is simply the 

Lebesgue space L'’(H). When to > 1 and r = 2, the space is a Hilbert space and is also denoted by 

For simplicity, denote by || • |j the L^(r2) norm, while other Sobolev norms shall be explicitly specified 
in subscripts, for example || • |l_f/i(o) and || • ~ note that in terms of norms there is no difference 

between and and hence the per i s om itted. For to < 0, denote by the dual space 

of Then the weak problem for Equation (1.5) can be written as: Find h, in a proper space to be 

specified later, such that for all t € (0, T) 


{dth, (j)) + a{h, (j)) = (C, (I)), for all </> e 
with the form a{h, 4>) = 7 (V/i, V0) + e^{Ah, Acj)) — ($(|V/i|)V/i, V0), 


( 2 . 2 ) 


where (•, •) denotes the duality pair, or the inner-product on H if both parties involved lie at least in 


Lemma 2.1 states that 0 < $(|V/i|) < C, thus the nonlinear term ($(|V/i|)V/i, V;/)) in (2.2) is well-defined as 
long as h and (j) are in 

Definition 2.6. We say h : H x [0, T] — >■ M is a weak solution to (1.5) if it satisfies 


1. he and dth € £^(0, T; 


2. Function h satisfies the weak formulation (2.2) almost everywhere for t G (0,T); 

3. h{-,0) = ho{-) almost everywhere in H. 

2.2.1. Semi-discrete Galerkin spectral approximation. Here we define the semi-discrete Galerkin 


spectral approximation to (2.2). For any given x S fl = (0, Li) x (0, L 2 ), denote x = 27r[a;i/Li, X 2 /L 2 y 


(0, 27r) X (0, 27r) . For a given positive integer N , define the index space Fj^j — G with 0 < 6,6 < 
N and ^ 7 ^ 0} and a discrete space on H by 

Hm = span{\, cos^ • x, sin| • x, for all | G Tat}- 

The space Fb^ is H-periodic. Note that the spanning set of also forms an orthogonal basis for Ha? under 
the inner-product. After proper ordering and normalizing, we get an orthonormal basis denoted by 

{fi, for i = 1,... ,S(77)}, where S(A^) = dimFlj^. Denote by P/v the projection onto Hat- We will seek 
the N-th. semi-discrete Galerkin spectral approximation to Equation (2.2) in the space as following: Find 
hj^r = h]S!^i(t)(j)i satisfying h]s[{-,0) = PNho{-) in D and 


{dthN,4>) = (C, </'), for all (j) G Hn, t G (0,T]. 


(2.3) 


We point out that sin ce 1 G H^!, the operator P^i maps mean value free functions to mean value free functions. 
By setting </> = 1 in (2.3), one has dt/iAf = 0. Gombining the above, we know that h^i, if existing, is mean 


value free for all t G [0, T] as long as C, and h^ are mean value free. 

Before establishing the well-posedness of the Galerkin spectral approximation (2.3), we first state two 
technique lemmas from m, with a little extra obvious facts. 

Lemma 2.7. For all ^ G FIp^,,{fl), one has 


2 

||V0f < m IIA0II and • 

i J = 1 


Moreover, if (j) is mean value free on by the Poincare inequality one has 

CU\\<\\Vf\\<^\\Acj^\\, 












where C is a positive general constant depending only on fl, and consequently ||^i>||_H'2(n) < 
Lemma 2.8. For any integer m > 0 and ( j ) G one has 


N—¥oo 


||P/V0||//m(Q) < \\(j)\\H’^{Q), and Pn4> - 4> strongly in 

Moreover, a direct calculation using Fourier series shows that 

U-PN^\\HHn)<CN-(^-^^U\\Hr^^n): for 0 < j < m, 

where C is a positive general constant. 

Next we prove the existence, uniqueness, and regularity of the Galerkin spectral approximation. 

Lemma 2.9. Assume that h^ G and C G ^^(0, T; then for each integer N > 1, there 

exists a unique semi-discrete Galerkin spectral approximation h^ satisfying pTM). The solution h^ has bound 


II^Af||L°°(0.T;L2(a)) + \\hN\\L^(0,T-H^(n)) < C'dl^oll, IICIlL2(o,T;ffp-r(n)))’ 


(2.4) 


where C'dl/ioH, ||C|li 2 (g -^.^-2 is a positive constant depending on ||/io|| ond ||CIIl 2 (q 2^.^-2 . Moreover, 

if ho G iLpg^(r2) and C G L^(0, T; L^(f2)), then we also have the following bound 

\\hN\\L’=°{0,T-H^{a)) + \\hN\\L'^{0,T-H‘i{a)) + ||3t/lAr||L2(o,T;L2(n)) < C'd|/lo||_f/2(Q), ||C||L2(0,T;L2(n)), (2.5) 

where C'(||/io||i/2(o), ||(||i2(o,r;L2(n)) is a positive constant depending on \\ho\\m(a) and ||CIU2(o,T;L2(n)• 

Proof. We follow the proof of Theorem 4.1 in [13], with some modifications on terms involving the ES 
current Jes- By setting cj) = (j)j ioi j = 1,... ,S(iV) in Equation (2.3) and using the orthogonality of basis 
functions, we get a system of ordinary differential equations 


dthN.jit) = /j(C(t),/iA.i(t), ■ ,/iA,H(A)(0)> for all j = 1,... ,^(fV). 


( 2 . 6 ) 


Condition /iAr(-, 0) = PNho{-) actually sets the initial condition h^jiG) = {ho,(j)j), j = 1,..., S(A^) for System 
(2.6). A standard procedure to prove global existence and uniqueness of the solution to (2.6) is to first get 
local existence and uniqueness by the Picard-Lindelof theorem, i.e., by showing that fj are locally Lipschitz, 
and then prove that the solution is bounded for t up to any given T/v < T. We first argue that all fj are locally 
Lipschitz. This indeed follows immediately from Corollary and the fact that composition, summation, and 
product of locally Lipschitz functions are also locally Lipschitz. 

Now by the Picard-Lindelof theorem. System (2.6) admits a unique local solution for t from 0 up to a Tat. 
By setting (j) = h^it) in Equation (2.3) and using lemmas 2.1 2.7 and the Young’s inequality, one has 




■ 7 ||V/iArf+ e 2 ||A/i^f = [ $(|V/iw|)|V/i ^|2 da; + (C, M 

Jn 


<C|i/r^f+ -||AdA,f+ C||C||^-2(^) 


diiA/. 


N\' 


where C is a positive general constant. Combining the HA/iAr|P terms, multiplying the inequality by 2e 
and integrating against t, then using the fact that < 1 for t G [0,T], we have for all r G [0,T/v] 


I^a(-, t)||^-I- 27|| VdAr||^2(Q .r;L2(n)) + ^^11 ^^A|li2(0_T-;L2(Q)) 


<llhjy(; 0 )ir+Ce^^^ll(ll 




< ll^olP + C'e 


,2CT 


(n))’ 


(2.7) 


where the last step follows from Lemma 2.8 Thus one has 


Z(N) 


= \\hN{-,T)\\^ < \\ho\\^ + Ce 


2CT 


2=1 




for all T G [0 , Tat], 


i.e., the solution to System (2.6) is bounded at Tat < T as long as /iq G and C, G L^(0, T; iJpg^(n)), 

Hence a unique extension of the local solution to [0,r], i.e. the global solution, exists. Moreover, Inequality 


(2.4) follows immediately from (2.7) and Lemma 2.7 
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To prove Inequality (|2.5|), we set (j) = dth^^t) in (2.3) and use Lemma 2.2 to get 




^l|V/i7vf+ + f G{VhM)dx 


= {cdth^)<l\\cr+l\\dthNr. 


Combining the ||9t/i7v|P terms and then integrating against t give 

2 


Jl|9*/iiv||i.(o..;L^(0)) + il|V/.^(-,r)f + ^||Ah^(,r)f < I||0)f + ^||0)f 


+ 


G'(VhAr(-, r)) dx 


G(V/lAr(-,0))dx 


2 IICIli2(0,r;L2(n))> 


( 2 . 8 ) 


for all T G [0,T]. Using Lemma 2.4 and|2.7[ we have for all r G [0,T] 


G(VhAr(-, t)) dx 


< [ \GiVhN{;T))\dx< [ {C + G\VhN{-,T)\)dx 

J Q J Q 

< G + G||Vh^(.,r)f < G + G||h^(.,r)f + ^||Ah^(.,r)f, 


(2.9) 


where G is a positive general constant. Combining ( |2.8[ )-(2.9) and applying Lemma 2.8 as well as (2.4) give 
Inequality (2.5) except for the ||/iAr||L 2 (Q bound. 


Finally, we shall estimate the \\hN\\L^(oT-H*(Q)) bound in Inequality (2.5). Denote by d any first-order 
spatial derivative. Setting (j) = —d^h^ in (|2.3|) and using integration by parts to get 


1 d 
2M' 


5/ijvf + 7l! Va/Gvf + e^ll - [ d [$(| W|)Vh^] • dx = -(C, 

Jn 


<c\\cr + l\\d^h^r. 


( 2 . 10 ) 


By Lemma [m especially noticing that $'(|VhAr|) < 0, we have 

$'(|Vfejv|) 

IV/lAfl 

Then, applying Lemma [2.7| to dh^ and and using the Young’s inequality give 


d [$(|V/i7v|)V/iAr] • WdhN = $(|V/iAr|)|V9/i7VP + " ' ^dhNf < G|V5/7vp. 


/ d [$(|V/iAr|)V/iAr] • VdhNdx < G||Va/iA,f < G||5/iA,f || A5/iArf 

Jn 

<G||ah^f+ ^l|Aah^f. 

Substitute the above inequality into ( |2.10 ), integrate against t, and use ( |2.4[ ), one gets 

II^Af||L2(0,T;ff3(Q)) < G(||ho||ij 2 (Q), ||CIU2(o,T;L2(n))- (2-11) 

Now similarly, set (j) = A'^h^ in ( |2.3[ ), integrate against t and apply the Young’s inequality. Again by 
Lemma ^ and Sobolev embedding, we have 




'0 Jn 


V • (4>(|V/ijv|)V/iAr) A^/iat dxdt 


—^(11/1011^2(0), ||CllL2(0,T;L2(n))) + 

<C/(li/io|| ^2(0), \\C\\hi0,T ;L2(n))) + ^ 

AG(|j/io||^2(o), ||Clli2(o,r ;L2(n))) + C 


10 Jn 

nT , 


$'(|V/iAr|) 


(V/lAr)*V2/lAr(V/lAr) 


|v/gvI 


4)(|V/lAr|)A/lAr 


A^/iat dx dt 


|V/iAr| IV^/ttvI + |A/rjv| 


\A^hi\f \ dxdt 


'0 JQ 
-T r 

II/lAT II ^1,4(0) II/lAT II ^2,4(0) + l|/'-Jv|lff2(o) 

^2 /-T 


dt+^ [ ||A^llAr|P dt 

4 Jo 


AG(||/io||^2(o)) ||CIIl2(o,T;L2(o))) l|/iAf|lL2(o_7’;//3(o))) + ^ / l|A hj^W dt. 


This, together with Lemma 2.7 and Inequality (2.11), completes the proof of the lemma. □ 
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2.2.2. Existence and uniqueness of the weak solution. Now we state the main existence, uniquess, 
and regularity theorem: 

Theorem 2.10. Let and ( € L^(0,T; L^(Q)), then System (1.5) has a unique weak solution 

h e L~(0,T;i?2^,(f2)) nL2(o,T;if4^,(n)) with dth G L^{0,T; L^{n)). 

Proof. We start from proving the existence of the weak solution. Using Equation ( |2.5[ ) and a compactness 
argument, it has been proved in Theorem 3.1 of m that there is a subsequence of Galerkin spectral ap¬ 
proximations hjv converging to a function h G L°°(0, T; with dth G L^(0, T; in the following 

sense: 

Af-foo 

^ h in L^{0,T; strongly. 

The rest of the proof of the existence also follows directly from the proof of Theorem 3.1 in m- Due to the 
different J es term, here we only need to re-prove Equation (3.14) in [12]. That is to prove for any (f G 
and r] G Cp, T], 


[ {rj{t)VPN^,mVhN{-,t)\)VhN{;t))dt- [ {rj{t)V<f,<^{\Vh{;t)\)Vh{-,t)) dt 
'o Jo 


N—^oo 


-G 0. 


( 2 . 12 ) 


Indeed, by Lemma 2.1 2.8 and 2.9 the left-hand side of Equation (2.121 satisfies 

LHS< [ (ry(t)[VPjv0-V<))],4>(|VhAr(.,t)|)Vhjv(-,t))dt 
Jo 


[<i>(ivh^(-,t)i) - $(|vh(.,t)i)]vh(-,t)) dt 


(r;(t)V<^,$(|V/i^(-,t)|)[Vh,v(-,t) - VM-,0]) dt 


< C'll^llL“(o.T)ll-PAf</^’—'/'llr/qn) / \\'^dE\\L^(n) dt 

Jo 

+ C'll^l|L“( 0 .T)||V())||i 4 (f 2 )||/ljv — h\\L2(0^T-,H^(n))\\'^h\\L2(0^T-,L^{Q)) 

+ C\\v\\L^io,T)\\^nmn) [ W^hN-^hWdt 

Jo 


N—>-oo 


-GO, 


V/iAr(-, t) — V/i(-, t) 


< C{T)\\Pm4> — 011^1(0) + C{T)\\hN — h\\L2{0,T;H^{n)) 
where in above we have used the fact that 

|$(|V/i^(-,f)|)-$(|VM-,i)|)| < fsup|$'(s)|) |V/ivr(-,t)|-|V/i(-,t)| <C 

\s>0 ) 

This completes the proof of the existence for the weak solution. 

To prove the uniqueness of the solution, let g and h be two solutions of System (1.5) with initial con¬ 
ditions go G iLpg^(D), /ig G iLpg^(D) and right-hand side functions g G L^(0, T; L^(D)), C, G L^(0, T; L^(D)), 
respectively. Define w = g — h. According to the previous proof of existence, one has w G L°°{0,T; Hp^^{Q)). 
Clearly for all (j) G iLpg^(D) and t G (0,r], w satisfies 

{dtw, cf) + 7 (Vu;, Vcf) + e(Au;, A<^) - ($(|Vg|)Vg - $(|V/i|)V/i, Vcf) = (g - C, <^). 


Set 4> = w{-,t) in the above equation. By the mean value theorem. Lemma 2.2 Lemma 2.3 Lemma 2.7 and 
the Young’s inequality, one gets 

^ Iwf + 7llVwf + e^llAwf = -{VpGiyg) - VfG(V/i), Vw) + (g - C, w) 


2 dt 


= —(V^G(m)Vu>, Vw) + {f] — 

<G||Vu;f+ h-CII \\w\\ 


<-||Au;f+ G||u;f+ h-Cf, 
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where m is a vector between Vg and Vh determined by the mean value theorem. Then, by the Gronwall’s 
inequality, we have 


ll'^llL“(o,r;L 2 (n)) + llwll< (7(1 + )(ll5o ^ o \\ l ^( q ) + \\ i 1 CllL2(-Q_7n.2^2(-Q))). 

The uniqueness of the weak solution follows immediately from the above estimate. 

Finally, we prove the regularity of the solution, i.e., h S L‘^{0,T-, Hp^^{n)). By Equation (2.5), is 
bounded in L^(0, T; i7pg^(n)). Hence by passing to a subsequence if necessary, we have converges to h 
weakly in L^( 0 , T; i.e., 

ll^llL2(o,T;H-i(n)) < limsup ||/iAr||i2(o_T;//'‘(n)) < <7- 

N^oo 

This completes the proof of the theorem. □ 

Remark 5. Similar to Theorem 3.3 in m , one can achieve higher order regularity if h^ and C are 
smoother. Also, existence of weak solution ean be proved with lower regularity requirement on ho and C, if 
one uses different spaces in the definition of the weak solution together with a refined compactness result, as 
discussed in m- Here we skip these details in order to quickly get a functioning well-posedness result that 
allows us to immediately start investigating numerical methods for the new model. 


2.3. Bounds of the solution and the roughness indicator. As shown in the proof of Theorem 2.10 


by using the weak convergence of h^ to h, one can pass the ||/iAr||L°°(o,T,Er 2 (j’)) < (7 upper bound to h. This 
upper bound, proved in Lemma |2.9[ is just a rough estimate. The purpose of this section is to derive more 
accurate upper bounds of the weak solution h in certain norms or semi-norms. To this end, we first state the 
following lemma: 

Lemma 2.11. For all (p G that is mean value free, one has 


where C is a positive general constant that does not depend on (p. 

Proof When Jes is set to defined in (1.2), by the definitions of $(-)j one has 


0< f mvp\)\VcP\^dx = aiiq-p) f 
Jq Jq 


nip+ \'^(p\){q+\'^<p\) 


dx < ai{q-p)\fl\, 


and when is set to J 2 ,£;s defined in (1.4), by Lemma 2.7 and Young’s inequality, one has 


0 < [ ‘^>{\\/p\)\\7(p\^ dx = a 2 [ 
Jq Jq 




Q q + {p + q)\'^(l>\/p 


dx < 


a2P 


|V 0 | dx < 


a2P 




(2.13) 


< C 


a2pVW\ > 

p + q 


P + q Jq 


P + q 


V0II 


This completes the proof of the lemma. □ 
Now we can prove the following estimate: 


Lemma 2.12. Let h be the weak solution to (1.5), then one has for all 0 < to < t < T, 


l\\hi;t)f+l f W^h^dr+U r ||A/.f dr < C{t-to) + h\h{;to)f + ^||CI| 

«/ tn ^ tn 




where C is the same constant as defined in Lemma 2.11, which depends only on ak, p, q, \il\, e and 7 . 


Proof. Setting p = h in (2.2) and using Lemma 2.11 give 


7||V/if+ e2||Adf = [ $(|Vh|)|Vh|2da:+(C,/i) 
Jq 


IQ 

<c+j\\Ahr 


rllCII 


Hy+{Q) 


-IIAhf. 
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Combine common terms and then integrate from tg to t, this completes the proof of the lemma. □ 

Remark 6. Lemma 2.12 states that ||/i|p, II dr, and || dr have at most linear growth rate 
with respect to t, starting from any point 0 < to <T. This is globally better than exponential growth. 

Remark 7. In the authors have studied the evolution of the surface roughness indicator for the finite 
ES barrier with slop selection case, i.e., the ES currents is Jfss.es- Here, by Lemma 2.1!^ we immediately 
have a same surface roughness evolution bound for the ES current 3k,ES with k = 1,2, since by definition one 
has |r2|a;^ = |l/i|p for t G [0,T]. JVhen C = 0, one gets a global bound for the growth of uj(t): 


uj{t) < ^yC{t — to) + ui^{to) for all 0 < to < t < T. 


(2.14) 


Next, we aim at deriving a local bound for the growth of uj{t) that is better than the global bound when 
the value of u}{t) is small. To this end, we first point out that the upper bound in Lemma 2.11 is not very 
sharp when |V(/)| is small. For small |V(/>|, one shall consider the following alternative bound: 

Lemma 2.13. For all G that is mean value free, one has 




(2.15) 


where C is a positive general constant that does not depend on fi. 

Proof. By lemmas |2.1[ |2.7| and the Young’s inequality, one immediately has 




Remark 8 . Comparing to Lemma 2.11. Lemma 2.1d\ is better when |1V^|| is small (or consequently, by 
the Poincare inequality, when ||0|| is small). 

Now we can derive the following local bound for ujft): 

Lemma 2.14. Let h be the weak solution to (1.5) and assume C = 0 - Then one has for all 0 < to ^ t < T, 

(2.16) 


uj{t) < oj‘^{to)e^d-to) ^ 


Proof. Setting = h in (2.2) and using Lemma 2.13 give 
1 d .... 


2 dt 


7 |!V/if+ e2||A/rf = f mvh\)\Vh\^ dx < C\\hf + U\Ahf, 
Jq 4 


which implies 


-I 

dt ' 


<cM^ 


one gets the local bound (2.16). □ 


Solve the ordinary differential equation and use = || 

Remark 9. Generally speaking, when t — to is large, the global bound (2.14), with growth rate 0 { y/t — to ), 
is much smaller than the local bound (2.16), which is an exponential growth. However, when uj{to) ~ 0 and 
t — to is small, the local bound (2.16) becomes smaller than the global bound (2.14 ). An illustration is given in 
Figure 2.1 The main reason that we derive the local bound is to explain that the local growth ofuj(t), when the 
value ofuj{t) is small, is indeed “flat” rather than the “abrupt” growth prescribed by (2.14 )■ Later, this pattern 
can be observed in numerical results. Other intermediate estimates between the local bound and the global bound 
can also be obtained easily through interpolation. For example, one can prove thatujit) < (C{t — to)+w3(to))i. 
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Fig. 2.1. Difference between the global bound ^^6 local bound \2.16\f , illustrated in an example with to = 0, 

aj(0) = 0.02 and C = 1. The local bound is better when t is small. More importantly, the local bound prescribes a “flat” start-up 
process in the beginning of the evolution, when u]{to) is small. 


Such intermediate estimates may he more accurate to describe the growth rate in some stages of the evolution. 
Just as we know, to strictly depict the growth rate in different stages is still one difficult task. 

Remark 10. When the ES current is taken to be 33 ,es, we do not have an existence and uniqueness 
theory due to the unboundedness ofJ^^ES o-t = 0- However, if there exists a weak solution, then the global 
and local hounds in ( 2 . 14 ) and {2.16} will apply, because one can easily show that for all cf S 




|v«.| 


(q-p) 


dx < aslfi] = C, 


dx < 


1- P Jq 4 


■cur, 


which are parallel to the results in lemmas 2.11 and 2.13 Thus inequalities {2.14) and {2.16} hold. 


2.4. Energy functional. Define an energy functional associated with Equation (1.5) as follows 


E{h) = f . 
in . 


GUh) + J\Vh{^ + -\Ah\^ 


dx. 


Then the differential equation (1.5) can be written as dth = { — where denotes the Frechet 

derivative of the energy functional with respect to h. 

Note that by definition G{Wh) can be negative, which implies that the energy functional E{h) can be 
negative too. However, by using the Poincare inequality, Lemma |2.7| and choosing f3 carefully, one immediately 
has the following lower bound of E{h): 

Lemma 2.15. The energy functional has lower bound 


E{h) > -C, 


where C is a positive constant independent of h. 

Another important observation is: 

Lemma 2.16. Let { € L°°(0,T, L^(fl)), the energy functional E{h) satisfies 

j^E{h)<l\\Cr-l\\dthr, for0<t<T. 


Note that E{h) is non-increasing when C = 0- 
Proof Setting (f = dth in (2.2) gives 


\\dur + j^E{h) = {c,dU) < IICII \\dth\\ < i|icf + luthf. 


This completes the proof of the lemma. □ 
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We shall discuss a little more about the special case when C = 0. In this case, it is not hard to see 
that there is a one-to-one correspondence between the critical points of E{-), i.e., h satisfying ^E{h) = 


^E{h + T(l))\T^o = 0, and steady state solutions of (1.5|. Moreover, local minimums of E{-) give stable steady 


state solutions of (1.51. We also point out that = 0 is obviously a steady state solution in this case. By 


Lemma [2.4[ E{h) becomes a convex functional of V/i when 7 is large enough. In this case, there exists a 
unique global minimizer, which is h = 0. Note here h = C is not considered because of the mean value free 
assumption. 

Remark 11. It is worth to point out that when the ES current is taken to be the energy functional 

is still well defined since G 3 G (7(11^) although it is not differentiable at the origin. Moreover, the lower bound 
in Lemma 1^.151 is also true because 


Gs{Vh) = -as ln(g -p+ |V/i|) > -G\Vh\ -G> - Gp, 


for any positive constant fi. 


3. A semi-implicit fully-discrete numerical scheme. In this section we develop a semi-implicit 


numerical scheme for approximating the weak solution of (1.51, using the technique of convex-splitting. For 


illustrative purpose, the spatial discretization uses the Galerkin spectral approximation with discrete space 
iJjv presented in Section though we point out that the scheme and analysis also apply to other Galerkin 
approximations. The numerical scheme is stated below. We first split the form o(-, •) defined in (2.2) into two 
parts: 


a+{h, 4>) = j{\7h, Vfi) + e^{Ah, A(j)) + {\7pG+{\/h), V(^), 
a_{h,fi) = {VFG_{Vh),Vfi), 

which satisfy a{h, fi) = a+{h, fi) + a-(h, fi). 

Given St = where M is & positive integer, and define fi = i St, ior i = t),..., M. Denote by /ijy S Hpi 
and Cf = C{-,ti), for i = 0,... ,M the numerical approximation and the deposition rate, respectively, at U. 
Then the semi-implicit time discretization can be written as: 

1. Set h% = PNho] 

2. For i = 1,..., M, compute h\ by 

l ^N Jt-n ^(j)) + a+{E^,(l)) = (3.1) 


for all (f> G Hn. 

By the definition of G+ and G_ in Gorollaryj^ it is clear that a+(-, •) is a symmetric and coercive bilinear 
form, while a_(-,-) is nonlinear. Hence the scheme is uniquely solvable at each time step. Next, we consider 
the energy stability of the scheme. 

For simplicity, denote 


g+ifi)= f G+(vfi) dx, g. {cf) = f G.(vfi) dx, 
Jq Jfi 


for ( j ) G and fc = 1, 2. 

Lemma 3.1. For G we have 


g+{(())-g+{il}) < [ VFG+{v4>)-v{(j)-tjj)dx, 

Ja 

I VfG_(VV') • V(()>-V')da:. 
Jq 


g_{fi)-g_{fi)< 
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Proof. By the mean value theorem and the fact that G+ is convex, one has 

Jq 


G+(V(^)-G+(VV')-VfG+(V(/>)-V(^-V') dx 
VfG+(V(/i - siV(^ - V')) - VFG+(V^i>) • V((/i - V') dx 
V|G+(V(/. - S2V((^ - V'))(-siV((/. - V')) • V((^ - V') da: 


< 0 , 


where 0 < Si < S 2 < 1 are constants determined by the mean value theorem. The proof for G_ is similar. □ 
Define the discrete energy at each time step i = 0, • • • , M by 

= E{h^) = + ^11 AhVf + ^ G(Vh5v) dx. 

we know that E^ > —C where G > 0 is a constant independent of More over, we have 


2.15 


By Lemma 

the following energy stability: 

Lemma 3.2. The scheme \3.1^ is unconditionally energy stable in the sense of 

E^ < G-i + |||Cf - - ^l|V(dV - - ^||A(hV - 


'-AT 


for all i = 1, ■ ■ ■ , M. Note that when Cf = 0, one has i?* < E' 




Proof Set (j) = in {3.1) and use Lemma 

we have 


3.1 


as 


well as the fact a(a — 6 ) = — 6 ^ + (a — &)^), 


'i— 11|2 


+E‘- E‘-‘ + ||| V(ftV - ft?')f + 1^11 A(ft? - ft?‘)f < (C‘, ft? - ft?') 


<fiicii“ + AL^ 


N 


This completes the proof of the lemma. □ 


Remark 12. The above energy stability indeed also implies H'^ stability. By lemmas 2.1 and 2.5. one has 

E^ > j\\Ah],r - /3||Vdjvf - G; 3 |D| > ^11 A/ftVf - Gp/3\\Ah],r - Gp\n\, 

where fd can be any positive constant and Cp is a constant from the Poincare inequality. Choose j3 such that 
2 

Gpfd = one has 


+ Cp\n\ > -IIA/ft^f > G-WE^Wl.^^y 

Thus the scheme is also stable. Note here we avoid using to control because 7 can 

be 0 . 

Finally, we study the error estimate of the numerical scheme. Denote the error at each time step ti, for 
z = 0 ,..., M by 

E = hix,ti) - 


where h is the weak solution to (1.51 and Ej^ is the numerical solution. The proof of the following theorem is 
quite standard and we thus postpone it to Appendix 
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Theorem 3.3. (Error estimate). Let Hq G and ( G L^{0,T; L^{Ll)). Assume the weak solution 

to 1 1.5) satisfy hu G L^(0, ht G L^(0, {LI)) with mi > 2, and h G L°°(0, i„_i; {Ll)) with 

m 2 > 2. Then 


\\hf\\ < C'e'='*-(<5t + 7V-™i 


Remark 13. Although we used the Fourier spectral Galerkin method in the spatial discretization, similar 
results hold for other Galerkin spatial discretizations. 

4. Numerical Results. An important feature of the semi-implicit fully discrete scheme (113 is that, it 
has a linear implicit part and hence one only needs to solve a linear problem in every time step. This is a great 
advantage comparing to the numerical schemes for models using or 3pss,ES as the ES current, which 

inevitably require a nonlinear implicit part for stability purpose [13 m]. However, one may wonder whether 
the new model, although easier to compute, can still correctly capture the evolution of surface morphology or 
not. In this section, we will first answer this question by comparing the numerical results from the new model 
with numerical results from other models reported in |13j . We will also test the new model on a larger set of 
examples to examine its performance. 

Set Ll = (0, 27r)^ and choose the Fourier spectral Galerkin approximation as the spatial discretization 
in Scheme We pick this spatial discretization because it can be easily and efficiently implemented in 

Matlab using the build-in Fast Fourier Transform (FFT) tool. In all numerical experiments, set the size of 
spatial discretization to be N'^ with N = 128, i.e., in the physical space h is evaluated on a 128 x 128 grid 
and consequently in the frequency space h is approximated by 128 x 128 Fourier modes, where h denoted the 
discrete Fourier transform of h. For simplicity, assume C = 0- 

Although is not continuous at V/i = 0, we have shown in Remark]^ that it still possesses several 

nice properties including the most important convex-concave splitting property. Thus we are also interested in 
testing k = S numerically and comparing it with fc = 1, 2. To distinguish between different ES currents 
^ 2 ,ES and here again we shall adopt the subscript fc = 1,2,3 throughout the rest of this section. In 

the implementation, one has to deal with the calculation of Ja.ES when |V/i| = 0. Here we adopt a makeshift 
solution by setting |V/i| = max{|V/i|, 10“^®}. 

Next we shall consider proper choice of the parameters, such as a, p, q, 7 and e, in Equation ( |1.5[ ). To 
make a realistic choice, let us first recall how the model was built from physical laws. According to [5], set 


ai 


FfEps 

2 


7 = CoEFf, 


Oi2 = C(1 


q-p 

p 


= Ff{L,si)\ 


a3 = ai{q-p), P=j —, 



b 

Les ’ 


where Ff is the deposition flux per unit time, Lps is the adatom attachment length when descending a step 
(ES effect), Lisi is the typical island separation length, b is the typical step height, and Cdf the strength 
of the downward tunneling current. We start from setting Lisi = 0.25, Lps = 0.05, b = 0.017, Cdf = 0 and 
Ff =2, which gives 


ai = 0.05, 02 = 0.25, 03 = 0.017, p = 0.068, q = 0.408, 7 = 0 , = 0.0078. 


Later we shall perturb parameters 7 and a little bit to investigate their effect on the surface evolution. Note 


that Equation (1.51 is linear in terms of t 
the time scale. 


therefore scaling a, 7 and together is equivalent to changing 


Thus we do not plan to test the numerical scheme for different values of a. Also, because 
of the small value of a we currently pick, the surface evolution with respect to time appear to be relatively 
slow. Hence we have found that setting the time step size St = 0.01 is adequate to resolve the rich details of 
the evolution. Though we point out that one may choose any other time step size and the numerical scheme 
will always be stable as proved in Section [3 However, St should be small enough in order to attain certain 
accuracy. One may also consider adaptive time-stepping strategies such as the one proposed in m- 


4 . 1 . Example 1. We start from the initial condition used in 


ho =0.1 (sin 3x sin 2y + sin 5x sin 5y). 
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Surface evolution with this initial condition using other models have been numerically studied in details in 
m- Interestingly, our numerical results show that the new model produces highly similar evolution patterns 
as those reported in m, despite the different ES current Jes used in these models. Note all these models are 
constructed based on the same physical phenomena, thus the numerical similarity indicates that they have 
each individually models the microscopic movement relatively correct. 



Roughness, for 0<t<300 





Fig. 4.1. Evolution of roughness uj and energy E for Example 1. Left: early evolution; Right: entire evolution. 
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Fig. 4.2. Evolution of surface pattern for Example 1. 


The evolution of surface roughness uj and energy E is reported in Figure [TT] It seems that a steady state 
solution has been reached at the end, as the roughness and energy curves appear to be flat. We shall point 
out that although the graph only shows evolution for 0 < t < 300, the actual computation is done for a much 























































































18 
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H1 semi-norm, for 0<:t<1.5 H1 semi-norm, for 0<t<300 




H2 semi-norm, for 0<t<1.5 H2 semi-norm, for 0<t<300 




Fig. 4.3. Evolution o/||V/i||^ and ||A/i||^ for Example 1. Left: early evolution; Right: entire evolution. 


longer time period, in order to ensure that the roughness and energy curves stay flat at the end. The same 
holds for all numerical results reported in this section. Contour plots of the solutions at different time steps 
in Figure 4.2 suggest that the case fc = 3 suffers from numerical round-off errors probably introduced by the 
crude treatment of |V/i| = 0. 

In all cases, notice that the energy in Figure |4.1| is non-increasing, i.e., the scheme is energy stable. 
Comparing with reports in |13) . they also share the following similarities: 


1. The roughness w drops in the beginning and then starts to increase, which has been described as a 
“rough-smooth-rough” pattern in [13]. 

2. The evolution goes through several “flat” stages, with each “flat” stage corresponding to a relatively 
stable surface pattern in the coarsening process. Similar “flat” stages have been reported in [T5] . 

3. We draw the surface image at different time, and report them in Figure |4.2[ Note that for either 
A: = 1 or fc = 2, the coarsening process evolves through three very different patterns, including the 
final steady state solution. We point out that these three patterns have exactly the same structure 
as the stages reported in [T3| . Also, there is no structural difference between the surface evolution for 
k = 1 and k = 2. The case fc = 3 is slightly different than fc = 1, 2 as the solution obviously is smeared 
by artificial round-off error in the middle of the evolution, which we suspect is introduced through 
the crude treatment of |Vh| = 0. Recall that the theoretical well-posedness and error estimate do not 
work for fc = 3. However, it is interesting to see that that numerical scheme for fc = 3 remains energy 
stable, and its solution eventually converges to a steady state of the same pattern as the solutions for 
fc = 1,2, with a small phase shift. In examples to be given later, we will see that there exist multiple 
types of steady state solutions for our model problem. But for a given initial condition, the behavior 
of solutions for fc = 1,2,3 seems to be similar and all three converge to the same type of steady state 
solution. 

4. We also point out that our numerical results have a longer evolution time length comparing to results 
in |l,3j . This is because we have picked small values for a and , which result in a slower coarsening 
process. 

The squares of semi-norms, i.e. ||V/i|p and ||Afc|p, of the solution for Example 1 are reported in Figure 
|4.3| The reason why we report squares of semi-norms instead of the semi-norms is that the squares instead 
of the semi-norms are components of E{h). One immediately notice that ||V/i|p grows in almost the same 
pattern as oj, while ||A/i|p appears to have sudden drops at the transition between “flat” stages. 

We also performed perturbation tests on parameters 7 and . Only the perturbation for fc = 1 is reported 
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Fig. 4.4. Effect of small perturbations of e^ and 7 on roughness for Example 1, with k = 1. 



Fig. 4.5. Effect of small perturbations of and 7 on surface patterns for Example 1, with k = 1. All solutions are drawn 
at time t = 250. 
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here since the test results are similar for /c = 2 and k = 3. The results are given in Figures 4.4 and 4.5 
which we draw the following conclusions: 


from 


1. Changing affects both the magnitude of ui and the time to reach the next “flat” stage. A smaller 
gives larger surface roughness at steady state, while slows down the evolution. This is reasonable 
as is the highest-order leading dissipation term in Equation (1.5). When is set to 0.512, the 

fourth-order dissipation term dominates the nonlinear ES effect, and the surface evolution behaves 
like a normal fourth-order dissipation, i.e., the roughness quickly drops to 0 and stays there, as shown 
Besides, from Figure 4.5 one can see that the magnitude of the steady state solution h 


2 . 


in Figure 4.4 


gradually drops to 0 as increases. 

Changing 7 affects the magnitude of oj and slightly affects the pace of the evolution. Again, this is 
reasonable as —jAh is the secondary dissipation term in Equation (1.5). Recall that in subsection 
|2.4[ we have drawn the conclusion that when 7 is large enough, the energy functional E is a convex 
functional of V/i and hence has a unique global minimizer h = 0. This has been observed when 
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we increase 7 to 0.64 in the perturbation test. In this case, the roughness drops quickly towards 0, 
which indicates that the stabilizing downward funneling current is dominant and the thin film growth 
becomes a simple dissipative process. Surface patterns of the steady state solution given in Figure [4^ 
also show that h is nearly 0 for 7 = 0.64. Again, note that the magnitude of the steady state solution 
h gradually drops to 0 as 7 increases. 

3. Another noticeable fact is that, the sign of the steady state solution can get reverted when perturbing 


as shown by comparing the first row of Figure 4.5 with the steady state solution for = 0.0078 


in Figure However, changing 7 seems to only affect the magnitude of the steady state solution. 


4.2. Example 2. In the second example, we pick an initial condition with high frequency: 

/iQ = 0.01 * (sin SOa: sin 20y + sin 50a:: sin 50y). 


Several stages of the surface evolution are reported in Figures [46l|4.8| for k = 1,2,3. All solutions converge 
to the same type of steady state solution that is different from the one for Examples 1, which indicates that 
the steady state solution is not unique. Again, the case fc = 3 suffers from numerical round-off errors. The 


evolution of roughness uj and energy E is reported in Figure 4.9 




Fig. 4.7. Evolution of surface pattern for Example 2, with k = 2. 


What is interesting about this example is the fast reduction of frequency (smoothing effect) in the beginning 
of the evolution. The initial condition has a high wave number that almost reaches the largest resolution of 
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k = 3 and t = 500 


k = 3 and t = 2000 
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k = 3 and t = 0.02 


k = 3 and t = 50 


Fig. 4.8. Evolution of surface pattern for Example 2, with k = 3. 
Roughness, for 0<t<20 Roughness, for 0<t<2000 
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Fig. 4.9. Evolution of roughness lj and energy E for Example 2. Left: early evolution; Right: entire evolution. 


a 128 X 128 computational grid. The surface plot in Figures 4.6||4.8 has to start from t = 0.02 oi t = 0.03, 
because the surface plot at t = 0.01 looks completely “black” due to its high frequency components. Within 
a few time steps, all these high oscillation parts are quickly smoothed out. This is obviously the effect of the 
fourth order dissipation term 

After the initial smoothing process, the magnitude of the solution undergoes a dramatic increase. For 
example, in Figure 4.6 the magnitude of the solution increases from 10“^ at t = 0.04 to 0.05 at t = 5 and 
eventually to 2 at t = 2000, which indicates a typical island-forming or coarsening process. Moreover, the 
wave number of the solution keeps dropping as t increases, until it reaches the steady state solution. 

Finally, we point out that although the case fc = 3 keeps suffering from round-off errors, it eventually 
converges to a same type of steady-state solution as the one for k = 1 or k = 2, again with a phase shift. 


4.3. Example 3. In the third example, we pick 


ho = sin 2 x cos 3y. 
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The purpose of this example is to test a relatively smooth initial data with large magnitude. 

Roughness, for 0<t<100 Roughness, for 0<t<3000 
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Fig. 4.10. Evolution of roughness u! and energy E for Example 3. 
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Fig. 4.11. Evolution of surface pattern for Example 3. 


Results of Example 3 are reported in Figures |4.10| and |4.11| Here we point out the sharp drop of both 
roughness and energy in the beginning of the evolution, as seen from Figure [4.10[ Indeed, comparing to other 
examples. Example 3 start from a large w, which is probably the reason of the sharp drop. This also agrees 
with the “rough-smooth-rough” pattern analyzed in |13j . 

Again, in Example 3, we found that the solution for fc = 1 and k = 2 converges to exactly the same 
steady-state pattern, with the magnitude for k = 2 larger then for fc = 1; while the solution k = 3 converges to 





































































































23 


a slightly different pattern with a phase shift but almost the same magnitude as for k = 1. This phenomenon 
has also been observed for examples 1 and 2. By examining the roughness and the energy history, one can 
also see that towards the steady state solution, the roughness and energy curves for A: = 1 and k = 3 tend to 
stay close while the curves for k = 2 are away from them. 

4.4. Example 4. In the fourth example, we pick 

ho = 0.01(sin3a;sin2?/ + cos 50a; cos lOOy). 


which is a combination of a low frequency part with a high frequency part. 
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Fig. 4.12. Evolution of roughness and energy for Example f. Left: early evolution; Right: entire evolution. 


k = 1 and I = 0.05 k = 1 and I = 70 k = 1 and I = 200 



Fig. 4.13. Evolution of surface pattern for Example f, with k = 1, 


Example 4 has the richest evolution process, i.e., the largest amount of “fiat” stages, among all examples 
presented in this paper. Moreover, it is the only example we have found so far such that there is a significant 
phase shift between the steady-state solutions for fc = 1 and k = 2. Results for Example 4 are reported in 
Figures |4.12||4.15[ 
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k = 2andt = 0.05 k = 2andt = 100 k = 2andt = 400 



Fig. 4.14. Evolution of surface pattern for Example f, with k = 2. 

k = 3andt = 0.05 k = 3 and I = 30 k = 3andt = 100 



Fig. 4.15. Evolution of surface pattern for Example f, with k = 3, 


4.5. Example 5. In the fifth and the last example, we pick a completely random initial condition with 
values in [—0.5, 0.5]. The initial condition was generated in Matlab using rand and then saved in a file, in 
order to make sure that all tests start from the same initial condition instead of another random generation. 

Evolution of Example 5 is reported in Figures 4.16[|4.19 An obvious “rough-smooth-rough” pattern is 
observed in the beginning of the evolution, as shown in Figure 4.16 Correspondingly, one can see how the 
“rough” random initial condition is smoothed out in Figures |T 

More interestingly, we found that for all k = 1,2,3, the solution goes through almost identical evolution 
stages, as shown in Figures |4 .1 7||4 .1 ^ Such a high similarity has not been observed in previous examples. 


Appendix A. Proof of Theorem 3.3 We first introduce a few notations. Note that a+(-, •) is a bilinear 


form and is coercive on This allows us to define an a_|_-projection from to iJjv by 


per \ 


a_|_(P^i;, (/)) = a_|_(u, (/)) ior all (j) G Hjy. 


It is standard to show that 

\\{I - for all tjt e with s = 0,1,2 and m > 2, 

where Ca is a positive constant independent of (p but might depend on e. 
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Roughness, for 0<t<1.5 
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Fig. 4.16. Evolution of roughness and energy for Example 5. Left: early evolution; Right: entire evolution. 
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Fig. 4.17. Evolution of surface pattern for Example 5, with k = 1. 


For simplicity, denote = h{-,ti) for i = 0,1,..., M. Define 


e* = P+h^ - h]^, p^ = h^- P+h\ 


then one has h* = e® + p®. We further denote 


0 ® = dfh^ - 


h® - h®-i 


St 


By Taylor expansion and the Schwarz inequality, one has 


= \\hui- 


' dt. 


Subtracting Equation (3.1) from Equation (2.2) gives 

{e\cl)) + -h^-\cl))+a+{h\cl)) = a_(/i®^\(/.) - a_{h\(l)) for all S 

































































Here we have used the fact that a+(-, •) is a bilinear form. The above equation can be further rewritten into 
^(e* - + a+{e\(l)) = -{9\<p) “ </>) + “ a-{h\4>). 

By setting (f) = 26te^, one gets 

Help - ||e-i|p + He* - + 26t(^ix + 7)HVe*f + e^||Ae*f) 

= - 26t{e\ e*) - 2(p* - /9*-p e*) + 26t(^a_{h]:^\e^) - a_{h\ e*)^ 

A/1+/2 + /3. 

By the property of 0*, one has for any constant Ci 

h<^5twr+c,6tMr<^ r \\hu{-M'^dt+c,6t\wr. 
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And for any constant C 2 , 

/2 = -2(r iI-P+)htdt,e^)<2\\ f' (/-P+Kdt||||e*|| 

Jti-i Jti-i 

UI - P+)htr dt^ ' lie*II 

<7^ r \\(.I-P^)h\\^dt + C 2 St\\er 

*-"2 iti_i 

f <2 rti 

< ||/i,||l,™,(^)dt + C2<5t||e*f. 

By Lemma 2.7, one has CpIjV^H < ||A^|| for all (j) G where Cp is the coefficient related to the 

Poincare inequality. Denote by Co = maxjy + 7, C|>e^}, which is to protect against the case when y + 7 = 0. 
Then one has 

Co||V<^f <(x + 7)||V0f+ e2||A^||2. 

By Corollary]^ There exists a positive constant Cq such that |V|:.G_(m)| < Cq for all m G Therefore 
I 3 = 25t{VpG.iyh^^^) - VpG_(VC), Ve*) < 25tCG\\V{h^^^ - /i*)|| ||Ve*|| 


< V^I|V(h^-'-h*)f+ <5tGo||Ve*|| 

< ^ {II Ve*-if + II + II V(h* - G-^)f ) + (5t ( (x + 7)llVe*f + e^|l Ae*|| 


z ||2 


Go 


Next, note that 


^l|Ve-f < 5|^||e-‘|| ||Ae-|| < ^||e-'||> + Al||Ae 


(A.l) 


and 


and 


5^IIVP-'||“ = 5^||V(/- (A.2) 

Co Oo Oo 


ZStCl 


^|lV(G-h*-i)f = 


Go 




G, 


0 Jti-x 


Wht dt|p < 


36eci 


G 


0 Jti-i 


W^htfdt. 


When 1 = 1, one shall replace the estimates in (A.l) and (A.2) by a combined term 

5^l|V(ft»„-A»)f = 5^^l|Ve»|7 

Choose Cl and C 2 to ensur e (G i + C 2 )St < Combine all the above, and sum up for i = 1,... ,n, use 
the definition of e° and Lemma 2.8 one has 


f ^^||Ve°f 
^0 


1 / 2 \ 
d|e"f+<5t^ (x + 7)l|Ve*f+ -||Ae*f <||e' 

+ Cdt^ ^ll^ttlli2(o.t„;L2(n)) + ll^tllz,2(o,t„;pi(a)) 

+ GA ||ht|||2(o_t^;pmi(f2)) + Gt„A ^^ll^lli°°fo.t„_,:g"»2(-ap +Gi5fy^ ||e' 

n —1 

< ||g2 


2=1 
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where C is a general constant that may depend on (7i, (^ 2 , Cp, Ca, C'g, 7, but not on St or N. Then 
the result follows from the GronwalPs inequality and the triangle inequality. 
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